MRS second order response

In this notebook we derive a correction for the spectral leak present in MRS band 3A detector images (see MIRI-TN-00002-KUL written by Dr. Bart Vandenbussche).

In [1]:
from IPython.display import HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
The raw code for this IPython notebook is by default hidden for easier reading.
To toggle on/off the raw code, click <a href="javascript:code_toggle()">here</a>.''')
Out[1]:
The raw code for this IPython notebook is by default hidden for easier reading. To toggle on/off the raw code, click here.

MIRI-TN-00002-KUL section 2: The MRS spectrometers are designed to cover the respective spectral band wavelength ranges in the first spectral order. A diffraction grating will yield maximum constructive interference for light of wavelength λ in the first order at an exit angle α. At the same exit angle α maximum constructive interference will also be seen for light of wavelength λ/2 (second order), λ/3 (third order), etc... Without a proper pre-selection of the wavelength bandpass entering the spectrometer, (parts of) these higher order spectra will be superimposed on the first order spectrum.

Such an occurence is expected to happen in MRS band 3A due to a spectral leakage from the dichroic transmission in band 1B. This is shown below.

In [2]:
# import modules
import funcs
import mrsobs

import numpy as np
import scipy.interpolate as scp_interpolate
from astropy.io import fits
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
plt.style.use('presentation')
%matplotlib notebook

import warnings
warnings.simplefilter('ignore')
In [3]:
# MRS spectral band of interest
band = "3A"
In [4]:
# Define paths to data
workDir   = '/Users/ioannisa/Desktop/python/miri_devel/'
cdpDir    = workDir+'cdp_data/'
d2cMapDir = workDir+'notebooks/distortionMaps/'
lvl2path  = workDir+'FM_data/LVL2/'
MRSWaveCalDir = workDir+"MRSWaveCal/" 
MTSDir        = workDir+'notebooks/pyMTSSim/OUTPUT/'
MrsFilterTransmDir   = MRSWaveCalDir+"MrsFilterTransmissions/"
In [5]:
# Read the measured transmission curves from the csv files
# zeroth colum is wavelength [micrometer]
# first column is room temperature transmission
# second column is 7K transmission
col = 2
filterWave= np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_1a.csv", delimiter=";")[:,0]
D1A = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_1a.csv", delimiter=";")[:,col]/100.
D1B = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_1b.csv", delimiter=";")[:,col]/100.
D1C = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_1c.csv", delimiter=";")[:,col]/100.
D2A = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_2a.csv", delimiter=";")[:,col]/100.
D2B = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_2b.csv", delimiter=";")[:,col]/100.
D2C = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_2c.csv", delimiter=";")[:,col]/100.
D3A = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_3a.csv", delimiter=";")[:,col]/100.
D3B = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_3b.csv", delimiter=";")[:,col]/100.
D3C = np.genfromtxt(MrsFilterTransmDir + "fm_dichroics_3c.csv", delimiter=";")[:,col]/100.
In [13]:
system_transmission = D1A*D2A*(1-D3A)

plt.figure(figsize=(8,6))
plt.semilogy(filterWave,system_transmission,zorder=0,label='system transmission')
plt.hlines(0.66,11.5,13,'lawngreen',label='100% order 1',zorder=1)
plt.hlines(0.0002,5.6,6.8,'r',label='0.02% order 2')
plt.hlines(0.0002,4,5,'orange',label='0.02% order 3')
plt.xlim(4,13)
plt.ylim(1e-6,1)
plt.xlabel(r'Wavelength [$\mu m$]')
plt.ylabel('System transmission')
plt.legend(loc='lower right')
plt.tight_layout()

Remarks:

  • The spectral leak of 2.8% at 6.1 micron will affect the MRS spectrum at 12.2 microns.

The contribution of the two spectral orders is disentagled as prescribed in MIRI-TN-00002-KUL. The result is shown below.

In [7]:
# Populate dictionaries with the wavelength, pixel size and validity maps 
# (corresponding to the detector plane)

d2cMaps   = funcs.load_obj('d2cMaps_band{}_tr10pc'.format(band),path=d2cMapDir)
sliceMap  = d2cMaps['sliceMap']
lambdaMap = d2cMaps['lambdaMap']
alphaMap  = d2cMaps['alphaMap']
betaMap   = d2cMaps['betaMap']
sizeMap   = funcs.get_pixel_spatial_area(band=band,d2cMaps=d2cMaps)

# create a "valid" map with all the pixels inside the slices of a particular band
channel = int(band[0])
# slice numbers in the slice map of the distortion CDP for this band
sliceInventory = np.unique(sliceMap)
slicesInBand = sliceInventory[np.where( (sliceInventory >= 100*channel ) & (sliceInventory <100*(channel+1)))]

validMap = np.zeros(sliceMap.shape)
for ss in slicesInBand:
    s = int(ss - 100*channel)
    # construct a list of y,x coordinates of detector pixels belonging to slices of this band
    pixels = np.where(sliceMap == ss)
    validMap[pixels] = 1
In [8]:
# Populate dictionaries with maps of L_sky
ip_Lsky = {}
for BBTemp in ["400K", "600K", "800K"]:  
    tabLSky = fits.open(MTSDir+ "MTSEquivalentLsky" + BBTemp + ".fits" )[1]
    ip_Lsky[BBTemp] = scp_interpolate.interp1d(tabLSky.data["wave"], tabLSky.data["L_sky"]*1000., kind='cubic')
In [9]:
EXTSOURCE_pixMap = {}
for BBTemp in ["400K", "600K", "800K"]:
    EXTSOURCE_SCI_pixMap,EXTSOURCE_BKG_pixMap = mrsobs.FM_MTS_BB_extended_source(lvl2path,band,bb_temp=BBTemp,corr1='',output='img')
    EXTSOURCE_pixMap[BBTemp] = EXTSOURCE_SCI_pixMap-EXTSOURCE_BKG_pixMap
In [10]:
R_n1,R_n2 = {},{}
for relation in ['800K_600K','800K_400K','600K_400K']:
    Resp_order1,Resp_order2 = [np.full((1024,1032),np.nan) for i in range(2)]
    sel = (d2cMaps['lambdaMap'] !=0)
    # solve 2x2 system of radiometric equations
    Teff1,Teff2 = relation.split('_')[0],relation.split('_')[1]
    #-- coefficients
    a1 = EXTSOURCE_pixMap[Teff1][sel]
    a2 = EXTSOURCE_pixMap[Teff2][sel]
    b1 = ip_Lsky[Teff1](d2cMaps['lambdaMap'][sel])*sizeMap[sel]
    b2 = ip_Lsky[Teff2](d2cMaps['lambdaMap'][sel])*sizeMap[sel]
    c1 = ip_Lsky[Teff1](d2cMaps['lambdaMap'][sel]/2.)*sizeMap[sel]
    c2 = ip_Lsky[Teff2](d2cMaps['lambdaMap'][sel]/2.)*sizeMap[sel]
    
    Resp_order1[sel] = (-a1/c1 + a2/c2) / (-b1/c1 + b2/c2)
    Resp_order2[sel] = (-a1/b1 + a2/b2) / (-c1/b1 + c2/b2)
    
    R_n1[relation] = Resp_order1
    R_n2[relation] = Resp_order2
In [11]:
fig,axs = plt.subplots(2,3,figsize=(12,8))
axs[0,0].imshow(R_n1['800K_600K'])
axs[0,0].set_title('Resp n=1 800K/600K')
axs[1,0].imshow(R_n2['800K_600K'])
axs[1,0].set_title('Resp n=2 800K/600K')
axs[0,1].imshow(R_n1['800K_400K'])
axs[0,1].set_title('Resp n=1 800K/400K')
axs[1,1].imshow(R_n2['800K_400K'])
axs[1,1].set_title('Resp n=2 800K/400K')
axs[0,2].imshow(R_n1['600K_400K'])
axs[0,2].set_title('Resp n=1 600K/400K')
axs[1,2].imshow(R_n2['600K_400K'])
axs[1,2].set_title('Resp n=2 600K/400K')
plt.tight_layout()
In [71]:
# use RSRF of band 3A as reference
fringe_file = cdpDir+'MIRI_FM_MIRIFULONG_34SHORT_FRINGE_06.02.00.fits'
photom_file = cdpDir+'MIRI_FM_MIRIFULONG_SHORT_PHOTOM_06.03.02.fits'

fringe_img     = fits.open(fringe_file)[1].data        # [unitless]
photom_img     = fits.open(photom_file)[1].data*0.55   # [DN/s * pixel/mJy]

# plot trace in pixel slice
ypos,xpos = funcs.detpixel_trace(band,d2cMaps,sliceID=d2cMaps['nslices']/2,alpha_pos=0.)

def tick_function(X):
    # for plotting wavelengths and wavenumbers on the same plot (two x-axes)
    V = X/2.
    return ["%.2f" % z for z in V]

fig = plt.figure(figsize=(12,6))
axs1 = fig.add_subplot(111)
axs2 = axs1.twiny()
axs1.plot(lambdaMap[ypos,xpos],(photom_img*fringe_img)[ypos,xpos],label='RSRF 3A mixed orders')
for relation in ['800K_600K','800K_400K','600K_400K']:
    Teff1,Teff2 = relation.split('_')[0],relation.split('_')[1]
    axs1.plot(d2cMaps['lambdaMap'][ypos,xpos],R_n1[relation][ypos,xpos],label='RSRF order 1 {}/{}'.format(Teff1,Teff2))
    axs1.plot(d2cMaps['lambdaMap'][ypos,xpos],R_n2[relation][ypos,xpos],label='RSRF order 2 {}/{}'.format(Teff1,Teff2))
axs1.hlines(0,11.4,13.5,'k',linestyle='dashed')
axs1.vlines(12.18,-0.5,3.25,'r',alpha=0.4,linestyle='dashed')
axs1.set_ylim(-0.2,2)
axs1.set_xlabel(r'Wavelength [$\mu m$] n=1')
axs1.set_ylabel(r'Response $[(DN/sec)\; /\; (mJy/pix)]$')
axs1.legend(loc='best',fontsize=10)
axs2.set_xlim(axs1.get_xlim())
tickmarks = np.array([11.5,12.,12.5,13.,13.5])
axs2.set_xticks(tickmarks)
axs2.set_xticklabels(tick_function(tickmarks))
axs2.set_xlabel(r'Wavelength [$\mu m$] n=2')
plt.tight_layout()
In [75]:
plt.figure(figsize=(11,8))
plt.imshow(R_n2['800K_400K'],cmap='jet',vmin=0.,vmax=0.11)
clb = plt.colorbar()
clb.set_label('RSRF order 2 $[(mJy/pix)/(DN/sec)]$')
plt.title("Second order response / band 3A")
plt.tight_layout()

Remarks:

  • The first and second spectral order responses determined from the different pairs of observations agree quite well, though not perfectly. We have found that a non-zero linear term in the non-linearity correction is the primal cause for this offset between responses.

Save first and second order response to CDP

In [23]:
# specify version to be saved
version = "7B.03.02"
In [24]:
# save responses to CDP
oldCdp = fits.open(cdpDir+"MIRI_FM_MIRIFUSHORT_SHORT_PHOTOM_06.03.02.fits")
dq_def = oldCdp[4]
In [25]:
outDir = cdpDir+"CDP7/"
In [26]:
def makeHeader(detector="MIRIFUSHORT", channel="12", band="SHORT"):
    hdu0 = fits.PrimaryHDU()
    hdu0.header["TELESCOP"]="JWST"
    hdu0.header["INSTRUME"]="MIRI"
    hdu0.header["MODELNAM"]=("FM", "Instrument model name")
    hdu0.header["DETECTOR"]=detector
    hdu0.header["DETSETNG"]="N/A"
    hdu0.header["READPATT"]="N/A"
    hdu0.header["SUBARRAY"]="GENERIC"
    hdu0.header["SUBSTRT1"]= 1                                            
    hdu0.header["SUBSIZE1"]= 1032                                              
    hdu0.header["SUBSTRT2"]= 1                                                
    hdu0.header["SUBSIZE2"]= 1024                                           
    hdu0.header["FASTAXIS"]= 1                                               
    hdu0.header["SLOWAXIS"]= 2
    hdu0.header["CHANNEL"] = channel
    hdu0.header["BAND"]    = band
    hdu0.header["FILENAME"]= "MIRI_FM_"+detector+"_"+channel+band+"_PHOTOM_{}.fits".format(version)
    hdu0.header["DATE"]="2018-12-20"
    hdu0.header["VERSION"] = version
    hdu0.header["USEAFTER"] ="2000-01-01T00:00:00"
    hdu0.header["AUTHOR"]  ="Bart Vandenbussche, Ioannis Argyriou"
    hdu0.header["ORIGIN"]  = "MIRI European Consortium"
    hdu0.header["EXP_TYPE"]= "MIR_MRS"
    hdu0.header["REFTYPE"] ="PHOTOM"
    hdu0.header["DESCRIP"] = 'CDP-7 MIRI MRS response'
    hdu0.header["PEDIGREE"] = 'GROUND'
    hdu0.header["PHOTMJSR"] = 42.5
    hdu0.header["PHOTUJA2"] = 1000.0
    hdu0.header.add_history("DOCUMENT: MIRI-TN-00003-KUL issue 1.2")
    hdu0.header.add_history("SOFTWARE: MIRICLE ")
    hdu0.header.add_history("DATA USED: RAL FM data obsId 11282 .. 11287")
    hdu0.header.add_history("DIFFERENCES: 06.03.00 Detector plane calibration (new format)")
    hdu0.header.add_history("DIFFERENCES: 06.03.01 Corrected error in PIXSIZ extension")
    hdu0.header.add_history("DIFFERENCES: 06.03.02 Fringe flat 06.02.00 applied")
    hdu0.header.add_history("DIFFERENCES: 7B.03.00 Renumbered to CDP7 Beta / fixed USEAFTER")
    hdu0.header.add_history("DIFFERENCES: 7B.03.01 Introduced 0.55 factor offset.")
    hdu0.header.add_history("DIFFERENCES: Used 10% transmission slice mask.")
    hdu0.header.add_history("DIFFERENCES: 7B.03.02 Separated the first and second order diffraction grating response.")
    return hdu0
In [ ]:
hdu0 = makeHeader(detector="MIRIFUSHORT", channel="12", band="SHORT")

resp = R_pixMap["800K"]["2A"].copy()
err = dR_pixMap["800K"]["2A"].copy()
size = sizeMap["2A"].copy()

sel = (validMap["1A"]==1)
resp[sel] = R_pixMap["800K"]["1A"][sel]
err[sel] = dR_pixMap["800K"]["1A"][sel]
size[sel] = sizeMap["1A"][sel]

sel = (( validMap["1A"] + validMap["2A"] )==1)
resp[sel] = resp[sel] # /ff["SHORTSHORT"][sel]

dq = np.full( validMap["1A"].shape, 2 )
sel = (validMap["1A"]==1)
dq[sel] = 0
sel = (validMap["2A"]==1)
dq[sel] = 0
dq = dq + BadPixelMap["800K"]["1A"]
dq = dq + BadPixelMap["800K"]["2A"]

hdu1 = fits.ImageHDU(data=resp, header=None, name="SCI")
hdu1.header["BUNIT"]= "DN sec^-1 mJy^-1 pixel"
hdu2 = fits.ImageHDU(data=err, header=None, name="ERR")
hdu2.header["BUNIT"]= "DN sec^-1 mJy^-1 pixel"
hdu3 = fits.ImageHDU(data=dq, header=None, name="DQ")
hdu4 = fits.ImageHDU(data=size, header=None, name="PIXSIZ")
hdu4.header["BUNIT"]= "arcsec^2"
hdulist = fits.HDUList([hdu0,hdu1, hdu2, hdu3, dq_def, hdu4])
hdulist.writeto(outDir + "MIRI_FM_MIRIFUSHORT_12SHORT_PHOTOM_7B.03.02.fits",overwrite=True)